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ABSTRACT 

New observational facilities are becoming increasingly capable of observing reflected light from 
transiting and directly imaged extrasolar planets. In this study, we provide an analytic framework to 
interpret observed phase curves, geometric albedos, and polarization of giant planet atmospheres. We 
compute the observables for non-conscrvativc Raylcigh scattering in homogeneous semi-infinite atmo- 
spheres using both scalar and vector formalisms. In addition, we compute phase curves and albedos 
for Lambertian, isotropic, and anisotropic scattering phase functions. We provide analytic expres- 
sions for geometric albedos and spherical albedos as a function of the scattering albedo for Rayleigh 
scattering in semi-infinite atmospheres. Given an observed geometric albedo our prescriptions can be 
used to estimate the underlying scattering albedo of the atmosphere, which in turn is indicative of the 
scattering and absorptive properties of the atmosphere. We also study the dependence of polarization 
in Rayleigh scattering atmospheres on the orbital parameters of the planet-star system, particularly 
on the orbital inclination. We show how the orbital inclination of non-transiting exoplanets can be 
constrained from their observed polarization parameters. We consolidate the formalism, solution tech- 
niques, and results from analytic models available in the literature, often scattered in various sources, 
and present a systematic procedure to compute albedos, phase curves, and polarization of reflected 
light. 

Subject headings: planetary systems — planets and satellites: general — planets and satellites: indi- 
vidual 



1. INTRODUCTION 

Observations of reflected light provide constraints on 
the scattering and absorption processes in planetary at- 
mospheres. In the solar system, geometric albedo spec- 
tra and phase curves have been used to infer chemi- 
cal compositions and cloud properties in several plane- 
tary and satellite atmospheres (Hansen & Hovenier 1974; 
Karkoschka 1994; Satoh et al. 2000; Sromovsky et al. 
2001; Irwin et al. 2002; Karkoschka & Tomasko 2011). 
Observational techniques are becoming increasingly ca- 
pable of detecting reflected light from exoplanetary at- 
mospheres, though currently available data are often 
limited to broadband visible photometry (Snellen et al. 
2009, Borucki et al. 2010, Demory et al. 2011). Our goal 
in the present work is to provide an analytic framework 
to aid in the interpretation of measurements of geometric 
albedos and phase curves of irradiated exoplanets. 

Several theoretical studies in the past decade have ad- 
vocated the use of reflected light observations for atmo- 
spheric characterization of extrasolar giant planets. Nu- 
merical models have been reported for directly imaged 
planets on wide orbits, as well as for closer-in giant plan- 
ets with very high incident stellar irradiation (Seager & 
Sasselov 1998; Marley et al. 1999; Seager et al. 2000; 
Goukenleuque et al. 1999; Sudarsky et al. 2000,2005; 
Stam et al. 2004; Sengupta & Maiti 2006; Burrows et 
al. 2008; Buenzh & Schmid 2009; Cahoy et al. 2011; de 
Kok et al. 2011; Kostogryz et al. 2011; Marley & Sen- 
gupta 2011). Atmospheres of extrasolar planets known 
to date, predominantly of extrasolar giant planets, expe- 
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rience a much wider range of equilibrium temperatures 
(Icq ~ 500 - 3000 K) than that of solar-system giant 
planets (Toq < 125 K). Sudarsky et al. (2000) suggested 
a classification scheme of albedo spectra for giant plan- 
ets, based on which highly irradiated planets (Toq ^ 900 
- 1500 K) are expected to have lower optical albedos due 
to strong Na and K absorption. At even higher tempera- 
tures, the presence of silicate condensates might increase 
the albedos. Albedo-based classification schemes for ex- 
oplanetary atmospheres were also discussed by Burrows 
et al. (2008), Cahoy et al. (2011), and Cowan & Agol 
(2010) as a function of different system parameters. 

Recently, observational constraints on the geometric 
albedos, visible phase curves, and potential polarization 
parameters, for several exoplanetary atmospheres have 
been reported. Some of the earliest constraints on geo- 
metric albedos {Ag) of extrasolar giant planets were re- 
ported for the non-transiting planets r Boo, with upper 
limits of 0.3 (Charbonneau et al. 1999), 0.22 (Collier- 
Cameron et al. 2000) and 0.39 (Leigh et al. 2003a), v 
And with an upper limit of 0.42 (Collier-Cameron et al. 
2002), and HD75289 with an upper limit of 0.12 (Leigh 
et al. 2003b). However, in recent years, with the dis- 
coveries of transiting exoplanets and with the advent of 
high-precision transit photometry, geometric albedo con- 
straints have been obtained for a substantial number of 
exoplanets. Rowe et al. (2008) observed the transit- 
ing hot Jupiter HD 209458b with the MOST satellite 
and reported a 3-sigma upper-limit of Ag < 0.17. Other 
hot Jupiters with Ag estimates include TrES-3 (Winn 
et al. 2008), CoRoT-lb (Snellen et al. 2009), CoRoT- 
2b (Snellen et al. 2010), HAT-P-7b (Christiansen et al. 
2010), Kepler-5b and 6b (Kipping & Bakos 2010; Desert 
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et al. 2011), HD 189733b (Berdyugina et al. 2011), and 
Kipping & Spiegel (2011), all of which have Ag uppcr- 
hmits less than 0.3. On the other hand. Kipping & Bakos 
(2010) and Dcniory et al. (2011) reported independent 
detections of a high geometric albedo for Kepler-7b of 
Ag = 0.38 ±0.12 and Ag = 0.32 ±0.03. Recently, several 
groups have also reported attempts to measure polar- 
ization parameters of exoplanetary atmospheres in the 
visible (Berdyugina et al. 2008, 2011a, b; Wiktorowicz 
2009). 

In this work, we derive an analytic framework with 
which to interpret observations of albedos, phase curves, 
and polarization of exoplanetary atmospheres. We gen- 
erally assume cloud-freqj homogeneous and semi-infinite 
atmospheres, and analytically solve the multiple scatter- 
ing problem for several different scattering phase func- 
tions. Detailed numerical models solve the general radia- 
tive transfer problem in a plane-parallel inhomogeneous 
atmosphere (Marley et al. 1999; Seager et al. 2000; 
Sudarsky et al. 2000; Burrows et al. 2008), with the 
assumption of radiative equilibrium for given chemistry 
and sources of absorption and scattering. Such models, 
however, involve significant computation time and con- 
vergence monitoring, and may be computationally pro- 
hibitive for formal fits to data, where several consecutive 
model evaluations are typically required. On the other 
hand, computationally efficient analytic models exist for 
several forms of scattering, if the atmosphere can be as- 
sumed to be homogeneous in the scattering albedcQ (uj). 
While the latter condition may not be exactly satisfied 
in giant planetary atmospheres across all temperatures, 
the analytic approach allows one to derive a represen- 
tative atmosphere-averaged scattering albedo from the 
observables, and has been successfully used in fitting ob- 
servations of solar-system planets (Kattawar & Adams 
1971; Sromovsky 2005; Schmid et al. 2006). 

Many researchers have derived analyticQ expressions 
for reflection coefficients, albedos, phase curves, and po- 
larization parameters for different forms of scattering 
in planetary atmospheres (Chandrasekhar 1950; Horak 
1950; Abhayankar & Fymat 1970; Kattawar & Adams 
1971; Van de Hulst 1981; Bhatia & Abhayankar 1982). 
Chandrasekhar (1950) and Horak (1950) reported ana- 
lytic solutions for reflected specific intensities from ho- 
mogeneous atmospheres for the cases of isotropic scat- 
tering, asymmetric scattering, and conservative Rayleigh 
scattering (w = 1.0), all for finite and semi-infinite at- 
mospheres. Chandrasekhar (1960) and Horak & Chan- 
drasekhar (1960) further generalized their solutions to 
address non-conservative Rayleigh scattering (w < 1) in 
semi-infinite atmospheres. Formulations using the scalar 
phase function for single-scattering do not take into ac- 
count the effect of polarization. Chandrasekhar (1960) 

^ We do co nsider the asymmetric scattering phase function (see 
Section I2.4.3|l which may be used as an approximation for par- 
ticulate scattering if the scattering asymmetry factor is known. 
However, such a treatment still does not address the full problem 
of Mie scattering by particulates. 

3 The scattering albedo is given by w = ascat/io-acat + crabs), 
where (Jscat is the single-scattering cross section and cTaba is the 
absorption cross section. 

By 'analytic' we mean solutions that can be expressed in closed 
form or which can be evaluated to arbitrary precision. This in- 
cludes results derived using ordinary (as opposed to partial) differ- 
ential equations. 



and Abhayankar & Fymat (1970,1971) reported a gen- 
eral solution including conservative and non-conservative 
Rayleigh scattering in homogeneous semi-infinite atmo- 
spheres using the full Rayleigh phase matrix, so that 
the Stokes parameters (Stokes 1852; Chandrasekhar 
1950,1960), and hence the polarization, could also be 
computed analytically. Subsequent studies have applied 
analytic approaches to a wide range of problems in plane- 
tary and exoplanetary atmospheres (Kattawar & Adams 
1974; Van de Hulst 1980; Bhatia & Abhayankar 1982; 
Sromovsky 2005; Schmid et al. 2006; Natraj et al. 2009; 
Kane & Gelino 2010). 

In this work, we present a systematic procedure to 
compute observables of reflected light for the different 
scattering phenomena under some basic assumptions. 
We consolidate the formalism, solution techniques, and 
results from analytic models that are available in the lit- 
erature, but are often in scattered and obscure forms. 
Building on the tradition of Chandrasekhar (1950,1960), 
we follow the H-function approach to represent emer- 
gent fluxes for different scattering phenomena. We con- 
sider cloud-free homogeneous semi-infinite atmospheres, 
scattering in accordance with different scattering phase 
functions. We consider both conservative (w = 1.0) and 
non-conservative (w < 1.0) scattering. We compute ge- 
ometric albedos and phase curves for Rayleigh scatter- 
ing, with scalar and vectorial phase functions, isotropic 
scattering, asymmetric scattering, and Lambert scatter- 
ing. We also compute polarization curves for the case of 
Rayleigh scattering. 

We provide a step-by-step procedure to obtain the 
disk-integrated emergent fluxes and polarization, with 
the appropriate angular transformations to the celes- 
tial reference frame. The phase curves and polariza- 
tion curves are also computed as a function of the mean 
anomaly, or a time coordinate, for given orbital parame- 
ters. We also provide analytic fits to the geometric and 
spherical albedos for Rayleigh scattering as a function of 
the scattering albedo. Thus, given an observed geomet- 
ric albedo at a certain wavelength, an average scattering 
albedo for the atmosphere can be estimated, which can 
then be used to predict other observables such as the 
phase curve and the Stokes parameters. 

In what follows, we first describe the formalism in Sec- 
tion [2] We present our results on phase curves and albe- 
dos for different scattering phase functions in Section [3l 
In Section |4l we study the dependence of polarization 
curves for Rayleigh scattering atmospheres on the orbital 
parameters. We summarize our results in Section [51 

2. METHODS 

Our goal in this study is to provide all the compo- 
nents of an analytic framework with which to interpret 
observations of reflected light from extrasolar planets. In 
this section, we briefly review some basic definitions of 
observables and describe the models used in this work. 

2.1. Definitions 

The spherical albedo (As) of a planet, which is the 
total fraction of incident light refiected by a sphere at all 
angles, for an incident flux irF, is given by: 

r i(a) 

As =2 :i^sina da, (1) 
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Fig. 1. — Planetary coordinate system. See Section 12.3.11 for 
description. 

where j(a) is the emergent flux and a is the phase an- 
gle, defined as the angle between the star, planet, and 
observer, whose vertex lies at the planet. As can also be 
expressed in terms of two other quantities of interest, the 
geometric albedo (Ag) and the phase integral (q), given 
by: 



ttF 



and (7 = 2 



jio) 



sin a da. 



(2) 



The quantity $(a) — j{a)/j{0) is the classical phase 
functior0. The planet-star flux ratio observed at Earth 
as a function of the phase angle is given by 



El 



^Ao 



1>(a), 



(3) 



where Rp is the planetary radius and a is the orbital 
separation. 

2.2. Orbital Phase 

Given a planetary phase curve as a function of the 
phase angle, the observed phase curve as a function of 
the orbital phase can be computed using Kepler's law. 
The phase angle (a) is related to the true anomaly {9) 
and the orbital inclination (i) by 

cos a = sin(6' + ujp) sin i, (4) 

where ujp is the argument of periastron. The true 
anomaly is related to the orbital phase via Kepler's law: 

(5) 



AI = E-esinE. 



^ The phase function, ^{a), should not be confuse d with the 
scatte ring phase function, p{cos ©) (sec e.g. Section 12.4.21 and 
[2A3t . 



Here, M is the mean anomaly, given by 27r(t — tp)/P; 
t is time; tp is time of pericenter passage; and P is the 
orbital period, e is the orbital eccentricity and E is the 
eccentric anomaly. E can be expressed in terms of the 
true anomaly as: 



sin£' 



sm 



eVT' 



1 + e cos 9 

The Kepler equation can be rewritten as: 



(6) 



t 



P 
2^ 



2 tan 




i9VT 



1 + e cos 9 



(7) 

Here, t = corresponds to pericenter passage. Thus, the 
planetary phase curve 4>(a) can be expressed as 

2.3. Emergent Intensity 

The emergent intensity, from a planetary atmo- 

sphere as a function of the phase angle determines all the 
reflected light observables, as is evident from Eqs. ^ & 
([3] ) . Theoretical computation of j (a) involves solving the 
radiative transfer equation while including the required 
sources of scattering and absorption in the atmosphere. 
Several numerical codes exist that solve the general ra- 
diative transfer problem of scattering in inhomogeneous 
atmospheres, i.e. where the scattering albedo (w) is a 
function of pressure (P) and temperature (T) in the at- 
mosphere, to compute albedo spectra and phase curves 
(Marley et al. 1999; Seager et al. 2000; Sudarsky et al. 
2000, 2005; Burrows et al. 2008). The sources of scat- 
tering typically include Rayleigh scattering, Mie scatter- 
ing due to condensates if present, along with a weaker 
contribution due to Raman scattering. For cloud-free at- 
mospheres, as are considered likely for highly-irradiated 
giant planets (Seager & Sasselov 1998; Sudarsky et al. 
2000), Rayleigh scattering dominates the scattering. In 
this work, we consider homogeneous (uniform w) and 
semi-infinite atmospheres for which analytic expressions 
can be derived for j{a) for several scattering conditions 
(Russell 1916; Chandrasekhar 1950; Horak 1950; Kat- 
tawar & Adams 1971), with particular emphasis on vec- 
tor Rayleigh scattering. 

2.3.1. Geometry 

The planetary coordinate system is shown in Fig[T] O 
is the center of the planet and OE and OS denote the 
directions from the planet to the Earth and the star, re- 
spectively. The angle EOS is the phase angle (a). CESG 
defines the planetary equator, and P and N are the poles. 
The longitude (C) and latitude (77) on the planet sphere, 
shown as dotted lines in Figlll are defined with respect 
to the great circles PEN and CESG, respectively. Q is 
an arbitrary point on the surface of the planet, shown 
for illustration, where a local coordinate system can be 
defined as shown. In this system, the z-axis is normal 
to the surface at Q. The incident ray from the star at 
the point Q is defined by the angles of incidence {9q, 0o) 
with respect to the local coordinate system. The x-axis is 
chosen such that 0o = 0, and the y-axis follows from the 
X and y axes assuming a right-handed Cartesian coordi- 
nate system. The direction of the reflected ray detected 
at Earth is defined by {9, cj)), as shown in FiglTJ The an- 
gles of the incident and reflected rays with respect to the 
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normal, a-nd 9, arc typically represented in terms of 
their direction cosines /iq = cos Oq and /i = cos 0. Vecto- 
rial quantities in the local frame at Q can be transformed 
to those in a global coordinate system centered at O by 
a rotation by an angle 7, as discussed in section 12.4.61 
below. 

2.3.2. Integration over the Visible Crescent 

At each phase angle in the orbit, in order to compute 
the total emergent intensity j{oi) from the planetary disk, 
we need to integrate the specific intensity over the illu- 
minated surface of the disk. In planetary coordinates, 
we denote the specific intensity of emergent radiation by 
7(77,^), where 77 and C are the planetary latitude and 
longitude, respectively, following the notation of Horak 
(1950) and Kattawar & Adams (1971). If /(t?, C) is known 
at each location on the disk, the disk-integrated emergent 
intensity at a given phase angle (a) is given by: 



df] sin" rj 



dC HvX) cosC. (8) 



a-7r/2 



The planetary coordinates (?/. (^) are related to the an- 
gles of incidence and reflection by 



^Q—smrj cos(C — a) 
1^1 = sim] cosC- 



(9) 
(10) 



The integral in Eq.® can be evaluated numerically 
by performing a coordinate transformation of Eq. ([S]) to 
obtain the following equation: 



+1 



j(a) = (£^£|±1) J ^ d^yr^ / d^ /(V,,e), 



+1 
1 



where, ip = cos 77, and 
2 



cos a + l 



cos a — 1 
cos a + 1 I ' 



(11) 
(12) 



where v = sin C- 

The double integral in Eq. pTj) can then be evaluated 
using standard quadrature methods as: 



J (a) 



(cos a + 1) 



^^WjWj/(V'j,^j), (13) 

i=i j=i 



where, Wi and Uj are the quadrature weights for the 
corresponding abscissae ipi and ^j, respectively (Horak 
1950; Kattawar & Adams 1971). In the present work, 
we use a 32-point Gaussian quadrature for the integral 
over ^, and a 32-point quadrature using Chebyshev poly- 
nomials for the ip integral. While the abscissae and 
weights for the Gaussian quadrature are available in stan- 
dard tables, the same for the Chebyshev polynomials can 
be evaluated simply by: ipi = cos[i7r/(ri + 1)], Wi = 
[TT/{n + 1)] sin^[i7r/(n + 1)]. 

2.4. Specific Intensity for Different Scattering 
Phenomena 

Before computing the total emergent intensity from the 
illuminated surface of the planetary disk, we need to be 
able to compute the emergent specific intensity at any 



arbitrary point on the planetary surface. Given an in- 
cident ray with coordinates (Oq, 0o) a-t a point on the 
planetary surface, as described in Section [2.3.11 our goal 
is to compute the reflected intensity in the direction (6, 
(/)). We shall refer to 9 and by their direction cosines 
jLt and HQ. 

2.4.1. Lambert Scattering 

One of the simplest models usually considered is that 
of an isotropically scattering Lambertian surface. The 
emergent specific intensity for such a surface (Chan- 
drasekhar 1960) is given by: 



(14) 



where (as defined earlier) uj is the scattering albedo and 
ttF is the incident stellar flux. As is evident from Eq. (fT4| . 
the reflected intensity is independent of the angle of in- 
cidence, and depends solely on the direction cosine (fi) 
of the observer. 

2.4.2. Isotropic Scattering 

Let us consider the case of a semi-infinite atmosphere 
scattering isotropically with a scattering phase function 
given by p(cos0) = w. For an incident flux {t:F) and 
a given angle of incidence (/ip = cos0o)i the emergent 
specific intensity at the top of the atmosphere (at optical 
depth T = 0) along a given direction (fj, = cos 6) is given 
by (Chandrasekhar 1950, 1960; Horak 1950)0: 



Hp) 



4 ^0 + ^ 



i/(M)i?(/io). 



(15) 



The function H{p,) satisfies the integral equation: 



1 



Hili) 



1 



■^{ti)dii 



(16) 

which can be solved by iteration until convergence is 
achieved. We found a reasonable initial condition to be 
i?(/i) = 1, and the integrals in the equation can be eval- 
uated numerically by standard quadrature methods; we 
used a 32-point Gaussian quadrature. vE'(/i) is known as 
the characteristic function, and is given by \I'(/i) = lo = 
constant, for the case of isotropic scattering. Thus, the 
emergent intensity at a given point on the planetary 
sphere can be computed for any angles of incidence and 
reflection. In order to compute the intensity detected by 
an observer at Earth, the emergent intensity in the di- 
rection of the observer is integrated over the illuminated 
crescent of the planet, as described in Section [2 .3. II 

2.4.3. Asymmetric Scattering 

The scattering phase function for asymmetric scat- 
tering is given by p(cos0) = a;(l -I- xcosQ), where Q 
is the angle between the incident and scattered rays, 
and X is the asymmetry factor, x € [0,1]. In analogy 

® We replace the notation 7(0; p.) of previous works (Chan- 
drasekhar 1950, 1960; Horak 1950; Abhyankar & Fymat 1970) with 
I{^^), since we are concerned only with the emergent intensity for 
which r = is implied. 
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with Eq. (|15p . the emergent specific intensity is given by 
(Chandrasckhar 1960): 



jF fio 



4 /io + A* 



V'(M)V'(Aio) - •^*(M)*(Aio) 



+ xJ{l - _ ^2)^(1) cos(0 - 



(17) 



where, = H{^){1 — cfi) and $(/^) = qfJ,H{^), and 

the characteristic functions with which to derive H{fi) 
and H^^\fi) are ^io[l + x{l - cj)^^] and jXuj{1 — yu^), 
respectively, q and c are given by q = 2(1 — w)/(2 — wao) 
and c = 2;a;(l — cj)ai/(2 — wao)j where ao and ai are 
the zeroth and first integral moments of H{fi): ao = 

/q H{^)d^ and ai ~ ^H{^)d^. 

2.4.4. Rayleigh Scattering - Scalar formalism 

Rayleigh scattering constitutes a dominant scatter- 
ing mechanism in planetary atmospheres, especially in 
the absence of Mie scattering by clouds. The scatter- 
ing phase function for Rayleigh scattering is given by 
p(cos0) = |a;(l-|-cos^ Q). This scalar form of the phase 
function does not account for the changes in polariza- 
tion in the reflected beam during multiple scatterings in 
the atmosphere. Nevertheless, it is customary in the lit- 
erature to adopt the scalar phase function for Rayleigh 
scattering for ease of computation (Marley et al. 1999; 
Sudarsky et al. 2000, 2005; Burrows et al. 2008). We 
discuss the emergent specific intensity with the scalar 
Rayleigh function in this section, and we discuss the full 
vectorial treatment in Section r2. 4. 51 

Horak & Chandrasckhar (1961) solved the case of 
a semi-infinite homogeneous atmosphere with a single- 
scattering scalar phase function given by 



p(cos9) = Lo + wiPi(cos9) + a;2P2(cos0), 



(18) 



where, the functions ^-"1(0;) and P2{x) are Legendre Poly- 
nomials of order 1 and 2, respectively, and w, Wi, and 
W2 are constants. The emergent specific intensity for 
the case of Rayleigh-like scattering in a semi-infinite at- 
mosphere can be obtained from their results by setting 
oji = and L02 = w/2 in Eq. to obtain the scalar 

phase function 



p(cos e) = -w(l + cos^ 9). 



(19) 



We note that the closed form expressions for the spe- 
cific intensity summarized in Section VIII of Horak & 
Chandrasckhar (1961) have some errors. Therefore, we 
choose to re-derive the specific intensity and obtained the 
following form for the solution, similar to Eqs. (|15p and 



4 Mo + 



$(Ai)$(Mo) + ^^(//)V'(mo) 



- |w'oy(r^M^)(l^^^'HM)i?^'HMo)cos((/.-0o) 



r(2) 



(20) 



where, $(m) = nH'^°\fi}^i{n) and V(m) = H'-°^{h)[3 + 
V'i(/i)]- H'-^^fi), H^^'^ifi) and iJ^^H/i) are H-functions 
with characteristic functions given by ^aj[3 — (2 — w)/i^-|- 
3(1 — a;)/x^], |a;(/z^ — A**), and ■^uj{1 — fj,'^)^ , respectively. 
$i(/i) and ipiifi) satisfy the coupled integral equations: 



'i>i(M) 



1 ,,2 



/iAi'$i(M)*i(/^') 



ff(o)(/x) 2 Jo + 

+ i(3 + Vi(A^))(3 + ^i(Ai')) dii' (21) 



and 



V'i(Ai) + 3 = 



3-Ai2 



i/(o)(Ai) 2 Jo + 



AiAi'4'i(Ai)$i(Ai') 



+ ^(3 + V^i(A*))(3 + Vi(Ai'))l (3 - Ai'')V- (22) 



The integral equations ([2T|) and ([22]) can be solved simul- 
taneously, for 3>i(At) and ipiifJ,), by iteration until conver- 
gence is achieved. We used 32-point Gaussian quadrature 
for the integrals. 



2.4.5. Rayleigh Scattering - Vector formalism 

In this section, we describe in detail the formalism 
with which to obtain the intensities of refiected radiation 
in directions perpendicular and parallel to the plane of 
incidence of a homogeneous and semi-infinite Rayleigh 
scattering atmosphere with arbitrary scattering albedo 
(oj). We follow the methodology of Abhyankar and Fy- 
mat (1970), who used the full Rayleigh phase matrix to 
develop an analytic formalism for scattering based on the 
H-funetion approach of Chandrasckhar (1950). 

The different components of reflected intensity, which 
are essentially elements of a modified Stokes vector 
(Stokes 1852; Chandrasckhar 1950,1960), are given by: 



Iiifi, 4>)=ujD [Al + Bi cos(0o - 0) + G cos[2((/)o - 0)] a^o^^ 
Iriix, (j))=UjD [Ar + Cr cos[2((/)o - <j))] ^0^^ (23) 
U in, (f)) = ujD [Bu sin(0o - (j)) + Cu sin[2(0o - 0)] ^o^^ 
and V{^l, 0) = . 



The Stokes vector is a four-element vector which de- 
scribes the intensity and polarization of a beam of light, 
and is represented as\ = [I Q U V]. Here. / is the total 
intensity, and Q and U are components of intensities with 
linear polarization, and V is the intensity with circular 
polarization. In terms of the modified Stokes parameters 
in Ea. (P5|) . I = Ii + Ir and Q ~ Ii — Ir, where // and 
Ir are the components of intensity parallel and perpen- 
dicular to the plane of incidence. The different terms in 
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Eg. ([25]) arc given by: 

Ar = iV3(/^)[A^l(M0) + + A^4(/x)[iV2(/io) + iV4(/Zo)] 



3 



(24) 



CjJ ~ 2^Cr 

and 

The functions Ni{fj,) and are computed along 

the hues of the general H-function approach of Chan- 
drasekhar (1950,1960). The A^j;(/i) satisfy the matrix 
equation: 



N{fi) = M{fi) + -ujnN(fi) / N^(^')M(m' 



d^l' 



(25) 



where. N 
and 

M = 



^1 N2 

N3 



N-^ is the transpose of N, 



" Ml 


M2' 


_^/3 - 


M3 


Mi 


2 



1 



\/2(l - 




(26) 



Eq. (|25p can be solved numerically by iteration until 
convergence is achieved. A reasonable initial condition 
is N = M. The integration on the right-hand-side of 
Eq. ipsj) is easily evaluated using a 32-point Gaussian 
quadrature. 

The H-functions H^^^ and ij'^^ are solutions of Eq. 
(1161) for characteristic functions: 



and 



M/W(/.) = -^(1-/.2)(1 + 2m2) 
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(27) 



(28) 



2.4.6. Transformation from local coordinates to celestial 
coordinates 



The intensities in Section 12.4.51 are defined with re- 
spect to a differential area element in the local meridional 
plane, i.e. the plane containing the incident and refiected 
rays at a local patch on the disk. These intensities need 
to be transformed into the global plane of reference of 
the planet as a whole with respect to the observer. The 
transformation is obtained by rotating the Stokes vector 
in the local plane through an angle 7 given by: 

cos7=:sin?7 sinC/sin0, or sin 7 = cos?]/ sin 0. (29) 

The rotation matrix for transforming I to the celestial 
frame is given by: 



10 

cos 27 sin 27 

— sin 27 cos 27 

1 



(30) 



The transformed intensities are given by I' = L • I, 
with the result that: 

/' = / 

Q' = Qcos27 + L/sin27 (31) 
[/' = -Qsin27 + [/cos 27 
V' = V . 

Since, /' = // + /; and Q' = ![ - /; by definition, the 
transformed ![ and I'^ are given by: 



I'l = {I' + Q')/2 and /; = (/' - g')/2. 



(32) 



These new intensities of the rotated Stokes vector are 
then used in Eq.((TT]) to compute the disk-integrated fiux 
in the planetary frame. 

The degree of polarization [P] is defined as P = 
-f U"^ / 1. For edge-on orbits, however, the disk- 
integrated Stokes U vanishes due to symmetry in the 
north-south direction, in which case P = Q/I- The an- 
gle of polarization is defined as, 6p ~ ^ i&n^^ [U / Q) . 
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Scalar Rayleigh 
Vector Rayleigh 




0.4 0.6 
Scattering albedo, a 

Fig. 2. — Geometric albedos as a function of scattering albedo 
for different scattering phase functions. The red and blue curves 
in the main panel show geometric albedos for Rayleigh scattering 
using the full phase matrix and using only the scalar phase func- 
tion, respectively. The inset shows the percent difference between 
the two curves. The green and orange curves in the main panel 
correspond to isotropic and Lambert scattering, respectively. 



3. RESULTS 

In this section, we present geometric albedos and phase 
curves, and, for the vector Rayleigh case, polarization 
curves, of reflected light for different scattering phenom- 
ena. As discussed in Section [2^ we consider the follow- 
ing cases: (a) Lambert scattering, (b) isotropic scatter- 
ing, (c) asymmetric scattering, (d) Rayleigh scattering 
with scalar phase function, and (e) Rayleigh scattering 
using the phase matrix, in which case the polarization 
can be computed. The geometric albedo {Ag) and phase 
curve, $(a), for each case are obtained by the following 
steps: 

1. For a given a, compute disk integrated emergent 
flux, j{a.), using Ea.([T3|: 
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Fig. 3. — Geometric albedos as a function of scattering albedo for 
the asymmetric scattering phase functions. For each ui, geometric 
albedos are shown for a range of asymmetry factors, x, between 
and 1. 



• Loop over summations: i = 1 — > 7ti. j = 1 — > n. 
in Eq. 

• Transform each ^j) to the corresponding 

• Compute intensity /(V'i, Cj) = Hp■^ for the 
given scattering case, using Eqs. ((T5|) . ([T7|) . (PH)) . 
or (Hi). 



• In the case of Rayleigh scattering with the vector 
formulation: 

— Compute the different intensities of the 
Stokes's vector in the planet frame [Irip-, 
Iiifi, (/)), U{fi, 0), V{fi, 0)] using Eqs. 
and 

— Transform the intensities to the celestial 
frame using Eqs. (|3l]) and ([32|) . 

— Use the transformed intensities in comput- 
ing different components of j{a). 

• Compute the final sum with the appropriate 
weights as shown in Ea. ((I3|) . 

2. Repeat Step 1 for q = — > tt and derive Ag, $(a), 
and polarization. 

• Compute Ag = j(0)/7rF, from Eq. © 

• Compute the phase curve, $(a) = j{a)/j{0). 

• For vector Rayleigh, compute polarization (see 
Section Em]). 

3. To obtain phase curves as a function of time (t), 
compute a for given t (see Section [2?2l ) 

3.1. Lambert Scattering 

The emergent flux for Lambert scattering is fully an- 
alytic (Russell 1916; Chandrasekhar 1950,1960) and is 
given by: 

2 p Sin(Q;) + (7r-a) cos(a) ^ 

j{a) = -ujttF[ J , (33) 

o 7r 



where uj is the scattering albedo and ttF is the incident 
stellar flux. Thus, the geometric albedo and phase func- 
tion of a Lambertian surface are given by = |w, and 
$(a) = [sin(a) + {n — a) cos(a)]/7r. For a perfectly re- 
flecting Lambertian surface (i.e. w = 1.0) the familiar 
Ag = 2/3 is obtained, and the spherical albedo is given 
hy As = 1. As shown in Fig. [2l the Lambert phase 
function leads to higher Ag for all uj, compared to all 
the other scattering mechanisms considered in this work. 
The phase curve for Lambert scattering is shown in Fig- 
ure □ 

3.2. Istrotropic and Asymmetric Scattering 

The cases of isotropic and asymmetric scattering cor- 
respond to scattering phase functions that are linear in 
UJ, as discussed in Sections 12.4.21 and 12 .4. 31 The isotropic 
phase function is a special case of the asymmetric phase 
function, with an asymmetry factor (x) of 0. The geo- 
metric albedo as a function of the scattering albedo for 
scattering in accordance with the isotropic and asym- 
metric phase functions are shown in Figs. [5] and [31 re- 
spectively. As shown in Fig. [51 isotropic scattering leads 
to lower Ag compared to that due to the Lambert and 
Rayleigh phase functions, for all uj. Figure [3 shows the 
phase curves for isotropic scattering for different uj. On 
the other hand, as shown in Fig. [3l asymmetric scatter- 
ing with a; > leads to higher values of Ag, compared to 
isotropic scattering, which for x I are comparable to 
those obtained from the Rayleigh phase function. 

3.3. Scalar and Vector Rayleigh .scattering 

Rayleigh scattering off gaseous atoms and molecules 
constitutes a dominant contribution to scattered light in 
cloud-free planetary atmospheres. Here, we compute the 
observables for Rayleigh scattering for both the scalar 
and vector cases, and we compare the results. Figure [H 
shows the geometric albedo (Ag) as a function of the 
single-scattering albedo for the two cases, and Figure [4] 
shows the corresponding phase curves, <i>(a). 

The geometric albedos for Rayleigh scattering using 
the scalar and vectorial treatments are not the same, as 
shown in Figure [^i. For conservative Rayleigh scatter- 
ing {uj = 1.0), i.e. a perfectly Rayleigh scattering atmo- 
sphere, the scalar phase function yields the commonly 
used value of Ag = 0.75. However, using the full phase 
matrix for the same case gives Ag = 0.7977, which is 
the more accurate value. The inset in Figure [2 shows 
the fractional difference between Ag in the two cases as 
a function of scattering albedo. Using the scalar phase 
function for computing Ag in the present case yields val- 
ues lower than the correct Ag by up to 9%, depending 
upon the scattering albedo. The differences between the 
two cases are lower for smaller uj, i.e. for greater absorp- 
tion. 

The phase curves for both vector and scalar Rayleigh 
scattering for different values of uj are shown in the top 
panel of Figure [H At large phase angles {a > 110°), 
the phase curves for different uj are similar. However, 
for a < 110°, phase curves for lower ui have steeper 
gradients, for both the scalar and vector cases. For a 
given UJ, the differences between the phase curves for 
Rayleigh scattering obtained using the scalar and vec- 
tor approaches are minimal. As shown in Figure [4l for 
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Fig. 4. — Phase curves and polarization for Rayleigh scattering. 
The different curves in each panel correspond to the different scat- 
tering albedos shown in the legend, fn the upper panel, the solid 
curves show the phase curves for Rayleigh scattering using the vec- 
torial Rayleigh phase matrix, whereas the dotted curves show those 
with the scalar phase function which does not incorporate polariza- 
tion. The lower panel shows the degree of polarization (P), using 

the Rayleigh phase matrix. P is defined as P = y^Q^~+~c7^//, 
where Q and U are the two Stokes parameters for linear polariza- 
tion and / is the total intensity. For all the curves shown here, the 
orbit is assumed to be edge on (i = 90°), in which case U = and, 
hence, effectively P = Q/I. See Section [4] for results with different 
inclinations. 

low scattering albedos (w < 0.5) the phase curves ob- 
tained by the two approaches are nearly identical. As 
with Ag, the differences grow with w, but additionally 
the differences are now also a function of the phase an- 
gle. Even so, the maximum difference between the curves 
for w = 1.0 is only ~5%. Such differences in phase curves 
are not observable given the current precision of obser- 
vations (Snellen et al. 2010; Demory et al. 2011). 

The degree of polarization of reflected light for a tran- 
siting planet with the Rayleigh phase matrix is shown 
in the lower panel of Fig. |4l The polarization curves, 
as a function of the phase angle, are shown for different 
uj. For all ui, the polarization peaks at a ^ 90°. How- 
ever, as is well known in the literature (Buenzli & Schmid 
2009), the degree of polarization (P) increases with de- 
creasing Lu. As shown in Fig. [SJ for highly absorptive 
atmospheres, say for w = 0.1, while the total intensity 
of reflected light is very low, the maximum polarization 
{Pmax) is high, approaching 100%. On the other hand, 
for highly scattering atmospheres, Pmax can be as low 
as 30%, for oj ~ 1.0. Since both the geometric albedo 
and the polarization curve for a given object at a given 
wavelength are related to the scattering albedo, a mea- 
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Fig. 5. — Maximum polarization due to Rayleigh scattering as 
a function of the scattering albedo. The degree of polarization is 
given by P = \/Q^hU^ / 1 and is a function of the orbital phase. 
The vertical axis shows the maximum attainable P in an orbit for a 
given scattering albedo, shown as the horizontal axis. For circular 
edge-on orbits the peak polarization in attained close to quadrature 
phase, o ~ 90° . 
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Fig. 6. — Analytic fits for geometric and spherical albedos as a 
function of scattering albedo for conservative and non-conservative 
vector Rayleigh scattering. The red and blue circles show the 
geometric and spherical albedos, respectively. The black curves 
through t he cir cles show the corresponding analytic fits, discussed 
in Section 13.41 The inset shows the accuracy of the fits for each 
case. 

surement of one can provide constraints on the other. 

3.4. Fits to Geometric and Spherical Albedos 

It is useful to quantify the dependence of the geometric 
and spherical albedos on the single-scattering albedo, so 
that given an observational albedo measurement one can 
constrain the underlying absorptive properties of the at- 
mosphere. As is evident from Eq. p3p . for Lambert scat- 
tering, Ag and Ag are proportional to lo. However, for 
other scattering phase functions, such as that of Rayleigh 
scattering, the dependence is non-linear. Even though 
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Fig. 7. — Comparison of phase curves for different scattering 
phase functions. The phase curves for Rayleigh scattering (both 
scalar and vector) and isotropic scattering are shown for several ui 
values between and 1; higher phase curves correspond to larger 
uj. For Lambert scattering, the phase curve is independent of u). 
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Fig. 8. — Comparison between phase curves due to Rayleigh scat- 
tering (with vector formulation) and Lambert scattering. The solid 
curves in the main panel show three phase curves for the cases de- 
scribed in the legend. The multi-colored dashed curves in the inset 
show the ratios between pairs of solid curves, e.g. the magenta- 
green dashed curve is the ratio between the magenta and green 
solid curves. 



(l-0.15.s)(l-5) 
(1 + 1.05s) ' 



(36) 



where s = \J\ — uj. As shown in Fig. IH the fit for Ag is 
accurate to within 1% for w < 0.99, and to within ^ 3% 
for ijj ^ 1. For u) ~ 1.0, a better fit is obtained for Ag 
without the term (0.95 -I- 0.08w) in the denominator in 
Eg . ((35l) . The fit for Ag is accurate to within 2% for all 

UJ. 

The above expressions can be used to retrieve the scat- 
tering albedo from an observed geometric albedo, under 
the assumption of a homogenous, semi-infinite Rayleigh 
scattering atmosphere. The w thus obtained is a repre- 
sentative average scattering albedo at a given wavelength 
of the observation. 

3.5. Comparison between different phase curves 

It is instructive to compare the phase curves resulting 
from different single-scattering phase functions. Figure[7] 
shows the phase curves for scattering due to Rayleigh 
scattering, isotropic scattering, and Lambert scattering. 
While the Lambert phase curve is independent of w, the 
isotropic and Rayleigh phase curves vary with ui. We find 
that the Lambertian and isotropic phase curves predict 
up to a factor of ^2 higher phase functions compared 
to the Rayleigh phase curves for all w for phase angles 
below ^120°, with the maximum differences occurring 
for a ^ 70° — 90°. Furthermore, for Rayleigh scattering, 
higher w leads to higher phase functions for most phase 
angles. As shown in Figure O an w = 1.0 results in up 
to a factor of ~1.5 higher phase function compared to a 
phase curve with w = 0.3, and up to ~1.3 higher phase 
function compared to a Lambert phase curve. Conse- 
quently, observations at phase angles corresponding to 
maximum differences in the phase curves can potentially 
be used to constrain scattering phenomena in planetary 
atmospheres. The magnitude of differences between the 
phase curves due to the different scattering phenomena 
are comparable to the those between phase curves of 
planets with different atmospheric and orbital charac- 
teristics. For example, Sudarsky et al. (2005) reported 
theoretical phase curves of extrasolar giant planets as a 
function of orbital distance, between which typical differ- 
ences are comparable to those seen in the phase curves 
in Fig. El and Fig. [51 



Ag and As are expected to increase monotonically with 
w, the exact functional form can be non-trivial, but can 
be obtained by fits to numerical results, van de Hulst 
(1974) suggested the following functional form for Ag 
for a semi-infinite atmosphere with a homogeneous cloud 
layer: 

(l~0.139s)(l-s) 

(1 + 1.170s) ' ^ ' 



Ag 



where, s = a/(1 — w)/(l — ijjg) and g is the scattering 
asymmetry factor given by < cos 9 > . 

Inspired by Eq. (p4|) . we obtain fits to Ag and Ag for 
non-conservative Rayleigh scattering in semi-infinite at- 
mospheres, using the vector formulation, as follows: 



Aa = 0.7977 



(1 - 0.23s)(l - s) 
(l + 0.72s)(0.95 + 0.08a;) 



(35) 



3.6. Scattering Albedos from Model Atmospheres 

As discussed in previous sections, the geometric albe- 
dos, phase curves, and polarization of scattered light de- 
pend strongly on the scattering albedo. The scattering 
albedo is given by, lo = a gcat I {cr gcat + crabs), where agcat 
is the single-scattering cross section and Uabs is the ab- 
sorption cross section. As such, the scattering albedo 
is a strong function of wavelength (A), and of the at- 
mospheric parameters such as temperature (T), pressure 
(P) , and composition, all of which govern agcat and dabs ■ 
Figure [9] shows scattering albedos in the visible wave- 
lengths, for a wide range of T, and for P of 0.1 and 1 bar 
which are representative of the pressure layers probed 
by visible/near-IR observations of exoplanetary atmo- 
spheres (Burrows et al. 2008; Sharp & Burrows 2007; 
Madhusudhan et al. 201 la, b). Here, we assume cloud- 
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Fig. 9. — Scattering albedos from model atmospheres. Scattering albedo spectra are shown for a representative range of pressures 
and temperatures observed in irradiated giant planet atmospheres. The models assume hydrogen-dominated atmospheres, with solar 
abundances, in chemical equilibrium. The scattering albedo is given hy id = (Jscat / {crscat + <Jabs), where ascat is the single-scattering 
cross section and cTaba is the absorption cross section. For wavelengths (A) blueward of ~ 0.55 fim, significant absorption due to Na, K, 
and molecular speci es ca use low scattering albedos. For A < 0.55 fim, Rayleigh scattering becomes prominent, leading to high scattering 
albedos. See Section 13.61 for discussion. 




Fig. 10. — Phase curves as a function of time from periastron pas- 
sage for different orbital parameters. The y-axis shows the phase 
function. The orbital parameters are shown in the panels: eccen- 
tricity (e), inclination (i), and argument of periastron {u}p). In each 
panel, the solid red and blue curves correspond to vector Rayleigh 
scattering with ui = 0.3 and oj = 1.0, respectively. The dotted red 
curve corresponds to Lambert scattering, which is independent of 



free, solar-abundance atmospheres in chemical equilib- 
rium, at the given P and T. 

As shown in Fig. [5J the scattering albedos are signifi- 
cant at short wavelengths (A < 0.55), and decline steeply 
to very low values (a; < 0.2) for A > 0.55. The increase 
in oj with decreasing A is due primarily to Rayleigh scat- 
tering. Furthermore, a; is a strong function of tempera- 
ture. Extremely irradiated atmospheres with T > 2000 
K provide substantial molecular absorption, leading to 
low scattering albedos across all A. The lower tempera- 
ture end of irradiated atmospheres, with representative 
temperatures of 1000 - 1500 K has the highest w at short 
wavelengths. Finally, w also depends on the pressure, 
and is higher for P = 0.1 bar than for P = 1 bar, and 
is generally higher for lower pressures. Consequently, an 
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Fig. 11. — Ph ase curves as a function of time from periastron 
passage. See Fig. llOl for description of axes. The orbital parameters 
for the curves are shown in each panel. The left and right panels 
in each row show the effect of the argument of periastron, for the 
same eccentricity and inclination. The top and bottom panels in 
each column show the effect of eccentricity, for the same argument 
of periastron and inclination. In each panel, the solid red and 
blue curves correspond to vector Rayleigh scattering with uj = 0.3 
and OJ = 1.0, respectively. The dotted red curve corresponds to 
Lambert scattering, which is independent of uj. 

average scattering albedo derived from data in a given 
spectral bandpass is indicative of the chemical and ther- 
mal properties of the atmosphere probed by the obser- 
vations. 

3.7. Time Dependence 

In previous sections, we expressed the phase curves for 
different scattering phenomena as functions of the phase 
angle (a). However, an observed phase curve, as a func- 
tion of time, also depends on the orbital parameters, as 
discussed in Section [521 Phase curves due to Rayleigh 
scattering and Lambert scattering as a function of time 
from periastron passage are shown in Figs. [T0lll21 for dif- 
ferent orbital parameters. For Rayleigh scattering, we 
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Fig. 12. — Phase curves as a function of tim e fr om periastron 
passage for different orbital parameters. See Fig. llOl for description 
of the colors. The orbital parameters are shown in the panels. 

consider two examples corresponding to = 1.0 (con- 
servative scattering) and lo = 0.3 (non-conservative scat- 
tering). The Lambert phase curve is independent of uj. 
In ah cases, the peak of the phase curve corresponds to 
a phase angle of 0°, i.e. full illumination at opposition, 
and the minimum corresponds to a phase angle of 180°. 
Here, we assume the longitude of the ascending node (ft) 
to be 90°. 

The differences between the phase curves due to the 
different scattering phase functions can be used to distin- 
guish between the underlying scattering phenomena, de- 
pending upon the orbital parameters. The phase curves 
for vector Rayleigh scattering generally have steeper gra- 
dients than those for Lambert scattering. However, the 
absolute differences between the curves vary with the 
orbital phase, and are a strong function of the orbital 
parameters. For example, for the case of e = 0, i ~ 
90°, Wp = 90° in Fig. [10] (top-left panel), the planet-star 
flux contrast for Lambert scattering with uj — 0.3 can be 
up to a factor of 2 higher than that for the correspond- 
ing Rayleigh scattering case. In this case, the maximum 
difference occurs at an orbital phase corresponding to (t- 
tp)/P ~ 0.2. It follows that observations made around 
this phase, in comparison with those made near the peak 
of the phase curve, could provide a constraint on the scat- 
tering mechanism. Thus, given the orbital parameters of 
a planet, similar times can be predicted for orbital phases 
which could provide the best diagnostics between the dif- 
ferent scattering phenomena. However, in some cases, 
such as in the top-right panel of Fig. [Til the differences 
between the different phase curves arc only marginal. 

4. POLARIZATION AS A FUNCTION OF ORBITAL 
PARAMETERS 



The polarization curves in Section 13.31 assumed the 
planetary orbit to be edge-on, namely that the inclina- 
tion angle (i) was 90°. Such an assumption is reasonable 
for planets in the solar system and for transiting extra- 
solar planets. In this section, we relax that assumption 
and study the dependence of the polarization on the or- 
bital parameters: inclination (i), eccentricity (e) and the 
argument of periastron (tOp). For the case of edge-on 
orbits, the Stokes parameter U is zero due to the latitu- 



dinal symmetry of the illuminated surface of the planet, 
whereas the Stokes parameter Q is non-zero due to asym- 
metry in the longitudinal direction. For an orbit that is 
not edge-on, however, both Q and U are non-zero, lead- 
ing to a non-zero polarization angle, X = 5 tan^^([/ /Q). 

The polarization parameters provide strong constraints 
on the inclination of the planetary orbit, as is well 
known in the literature on binary stars (e.g. Rudy and 
Kemp 1978; Schmid 1992; Harries and Howarth 1996). 
Recently, measurements of polarization parameters are 
also being attempted for giant exoplanetary atmospheres 
(Berdyugina 2008, 2011a,b; but cf Wiktorowicz 2009). 
Here, we explore the dependence of the polarization pa- 
rameters on the orbital properties of a planetary system 
using our model formalism developed in previous sec- 
tions. 

The time-dependent Stokes parameters in the ob- 
server's frame can be derived using the following pre- 
scription (also see Schmid 1992): 

1. Given a mean anomaly (M), the true anomaly (0) 
can be obtained by solving the Kepler's equation, 
Eq. ©. 

2. The effective phase angle corresponding to a true 
anomaly of 9 in an inclined orbit is given by 



a = sin ^ Y sin^ tj cos^ i + cos^ ?y ; < < (37) 
and 



a = TT — sin Y sin rj cos^ i + cos^ 77 ; tt < 77 < 27r (38) 

where, rj ~ (6 + ujp) is the true longitude, ujp is the 
argument of periastron, and i is the orbital incli- 
nation. 

3. Given the a, the Stokes vector (S) in the plane 
of the orbit is computed using the prescription in 
Section [3l 

4. Finally, the Stokes vector (S) must be rotated into 
the frame of reference of the observer. The rotation 
angle is given by: 



7 = tan [cot(77) / cos(i)] 



(39) 



Here, we have assumed the longitude of the ascending 
node (fi) to be tt/2. The observed Stokes vector (S') 
is obtained by S' = L • S, where L is the transformation 
matrix given by Eq. ([30| . 

Model polarization curves as a function of the orbital 
inclination for circular orbits are shown in Figs.[T3lfc[T4l 
Here, we assume conservative Rayleigh scattering with 
the vector phase matrix discussed in Section 13.31 The 
Stokes parameters (L Q, and U), the degree of polar- 
ization (P = y^{Q'^ + U^)/I), and the polarization an- 
gls ix = ^ tan^^(J7/(5)) are shown as a function of the 
mean anomaly (M) for different inclinations. For edge- 
on orbits {i = 90°), the Stokes U parameter vanishes at 
all orbital phases, due to symmetry in the illuminated 
planetary area in the north-south direction, whereas the 
Stokes Q is non-zero. For all other inclinations (i < 90°), 
both Q and U are non-zero. On the other hand, the 
planetary phase curve (i.e. Stokes I as a function of M) 
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Fig. 13. — Polarization curves for conservative Rayleigh scattering (ui = 1.0) assuming circular orbits (e = 0). Various quantities involving 
the Stokes parameters (I, Q, U), integrated over the illuminated planetary surface, are shown as a function of the mean anomaly (M). 
The Stokes parameters are expressed as a fraction of the incident intensity. I is the total intensity, and Q and U are the two polarization 
parameters for linear polarization. The various quantities are shown for different orbital inclinations. For edge-on orbits (i = 90°), the 
maximum value of I in the orbit gives the geometry albedo, and U = due to symmetry in the north-south direction. The degree of 
polarization, given by P = \/Q^ + IJ^ / 1 , and the polarization angle x = 0.5tan~^(C//Q) are also shown as a function of the mean anomaly. 
For the i = 90° case, the irregularities in the degree of polarization close to M = 180° are unphysical, caused numerically by division by 
vanishingly small values of I. 



is steepest for edge-on orbits and gradually flattens out 
for smaller inclinations, leading to a uniform phase curve 
for a face-on orbit (i ~ 0°). Secondly, while the maxima 
of Q/I and U/I in the orbit increase with decreasing in- 
clination, the minima remain relatively unchanged. 
The variation of the degree of polarization (P = 

\/ {Q^ + U^)/I) SI'S a function of M for different orbital 
inclinations is shown in the lower-right panel of Fig. 1131 
A face-on orbit yields a constant polarization at all or- 
bital phases, and an edge-on orbit causes the maximum 
peak-to-trough variation with orbital phase. Strategic 
monitoring of P over the orbit provides useful diagnos- 
tics on the orbital inclination of the system. The peak 
polarization (Pmax) occurs at the two quadrature phases, 
mean anomalies (Af) of 90° and 270°, independent of the 
inclination. In addition, Pmax is identical for all inclina- 
tions. On the other hand, the minimum polarization 
(^min) is a strong function of the inclination, and ranges 
0° and for i = 90°. For circu- 



lar orbits, Pinin occurs at the conjunctions, i.e. orbital 
phases of and 180°. Consequently, the inclination of 
a planetary orbit can be constrained by monitoring the 
extrema in the polarization phase curve of the planet. 

The influence of the orbital inclination on the polariza- 
tion parameters is also apparent in the trajectory of the 
planet in Q - U space. Fig. [14] shows the planetary or- 
bits in the Q/I - U/I space for different inclinations. For 
each inclination, the trajectory is a double loop structure 
which becomes increasingly eccentric with increasing in- 
clination. For face-on orbits a double-circular trajectory 
is obtained, whereas for perfectly edge-on orbits a hori- 
zontal trajectory (with U = 0) is obtained. Consequently, 
the ellipticity of a planet's path in the Q/I - U/I space 
can be used to infer the orbital inclination of the system. 
The arrows in each panel show the variation of the po- 
larization angle with the mean anomaly, M. The colored 
circles show the position of the planet in Q/I - U/I space, 
and hence its polarization angle (x — 5 tan~-^{U/Q)), at 
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Fig. 14. — Polarization curves for conservative Rayleigh scattering {ui = 1.0) assuming circular orbits (e = 0). The trajectory of the 
planet in the Q/l - U/I plane is shown for different inclinations. The arrows show the position of the planet in the Q/l - U/I plane at 
different mean anomalies (M), spaced apart by 15°. Sample locations for M of 0°, 90°, 180°, and 270°, are shown by the colored circles 
described in the legend. See section [4] for discussion. 



four different times (denoted by the mean anomaly, M) 
in the orbit. The arrows show the same at 15° intervals 
between M of 0° and 360°. 

The polarization parameters for eccentric orbits for dif- 
ferent orbital inclinations are shown in Figs. [15] & [1^1 For 
purposes of illustration, we consider an orbital eccentric- 
ity (e) of 0.5, and an argument of periastron (wp) of 60°. 
Again, we set the longitude of the ascending node (f2) 
to 90°. The variation of the total intensity (/) with the 
mean anomaly (M) is governed by the Kepler equation, 
Eq.[5l which relates M to the phase angle (a), depending 
on e and to. Contrary to the circular case, the variation 
in the polarization parameters is asymmetric in time for 
eccentric orbits. However, similar to the case of circular 
orbits, the mean anomalies corresponding to the maxima 



and minima remain largely unchanged with the orbital 
inclination. Also similar to the circular case, the peak 
polarization in eccentric orbits is independent of the or- 
bital inclination (i), whereas the minimum polarization 
is a strong function of i. 

The trajectory of an eccentric orbit in the Q/I - U/I 
space is shown in Fig. [121 As in the circular case, the or- 
bit traces a double loop structure which becomes increas- 
ingly eccentric with increasing inclination. However, for 
the eccentric case the polarization angles are naturally 
shifted in time, as shown by the positions of the planet 
in the Q/I - U/I plane at the different mean anoma- 
lies. Consequently, by observing a polarimetric orbit of 
a planet in the Q - U space, the inclination as well as 
the eccentricity and the argument of periastron of the 



14 



o 




0.4 - 



-0.2 - 




100 200 300 

Mean Anomaly, M 



100 200 300 

Mean Anomaly, M 




Cl 



0.2 - 








300 



M 



100 200 300 100 200 

Mean Anomaly, M Mean Anomaly, 

Fig. 15. — Polarization curves for conservative Rayleigh scattering for eccentric orbits. See description in Fig. 1131 For illustration, an 
eccentricity (e) of 0.5 and an argument of periastron {iiip) of 60° are assumed. The longitude of ascending node (f2) is assumed to be 90°. 
For this eccentric case, the Stokes parameters vary asymmetrically with the mean anomaly, contrary to the circular case. However, for 
edge-on orbits the Stokes U still vanishes. 



planetary orbit can, in principle, be determined. 

5. DISCUSSION AND SUMMARY 

Observations are becoming increasingly capable of de- 
tecting geometric albedos, phase curves, and polarization 
from cxtrasolar planets in reflected light. We have pro- 
vided an integrated analytic framework to interpret such 
observations. Our approach makes accessible a compu- 
tationally efficient procedure with which to estimate an 
average scattering albedo (w) from an observable in re- 
flected light, such as a geometric albedo or a phase curve. 
We follow the H-function approach of Chandrasekhar 
(1950,1960) which has been used in subsequent studies to 
derive analytic solutions for scattered emergent flux for 
various scattering phenomena. We consolidate the for- 
malism, solution techniques, and results from analytic 
models available in the literature, but often scattered in 
various papers, and present a systematic procedure to 
compute albedos, phase curves, and polarization of re- 
flected light. 

In this work, we considered cloud-free, homogeneous, 
semi-infinite atmospheres scattering in accordance with 
different scattering phase functions. We consider both 



conservative {uj ~ 1.0) and non-conservative {uj < 1.0) 
scattering. We compute geometric albedos and phase 
curves for Rayleigh scattering, with scalar and vector 
phase functions, isotropic scattering, asymmetric scat- 
tering, and Lambert scattering. We also compute polar- 
ization curves for the case of vector Rayleigh scattering, 
and provide a step-by-step procedure to obtain the disk- 
integrated emergent flux and polarization, for a given 
scattering albedo and phase angle, with the appropriate 
angular transformations to the celestial reference frame. 
The phase curves and polarization curves are also com- 
puted as a function of the mean anomaly, or a time co- 
ordinate, for given orbital parameters. 

Rayleigh scattering is a dominant scattering mecha- 
nism in cloud-free planetary atmospheres. An accurate 
calculation of the geometric albedo (Ag) for Rayleigh 
scattering must include the effect of polarization, us- 
ing the full Rayleigh phase matrix for single-scattering, 
contrary to the customarily adopted scalar phase func- 
tion. For, conservative Rayleigh scattering (w = 1.0), 
Ag = 0.7977, as opposed to the familiar value of 0.75 
obtained using the scalar treatment. The scalar formu- 
lation leads to geometric albedos that are lower than the 
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Fig. 16. — Polarization curves for conservative Rayleigh scattering for eccentric orbits. The trajectories in the Q/I — U/I plane for 
eccentric orbits arc rotated with respect to those for circular orbits for the corresponding inclinations showrn in Fig. 1141 The arrows mark 
the mean anomalies (M) in the orbit, spaced apart by 15°. The locations for M of 0°, 90°, 180°, and 270°, are shown by the colored circles 
described in the legend. The orbital parameters for each case are shown in the panels. 



actual values by up to 9%, depending on oj. We com- 
pute albedos and phase curves for Rayleigh scattering 
for a range of w, and report analytic fits to the geometric 
and spherical albedos as a function of uj. Based on the 
polarization curves, we also compute the peak polariza- 
tion due to Rayleigh scattering as a function of oj. The 
emergent intensity and degree of polarization of reflected 
light for Rayleigh scattering depend strongly on u). The 
emergent intensity at full phase, and, hence, the geomet- 
ric albedo, both increase with cj. On the other hand, the 
degree of polarization (P) decreases with increasing w. 

Motivated by recent efforts to measure polarization of 
reflected light from exoplanets, we compute polarization 
curves for homogeneous Rayleigh scattering atmospheres 
over a wide range of orbital parameters using the vector 



Rayleigh phase matrix. We explore the dependence of 
the Stokes parameters on the orbital inclination and ec- 
centricity. We demonstrate how the inclination of a plan- 
etary orbit can be constrained by monitoring the extrema 
in the polarization phase curve of the planet, as well as 
the ellipticity of the orbit in the Q-U plane. 
We summarize the results of our work as follows: 

• We present an analytic framework to interpret ob- 
servables of reflected light (phase curves, geometric 
albedos, and polarization parameters) from extra- 
solar planets. Our analytic models assume cloud- 
free homogeneous and semi-infinite atmospheres, 
and we consider several different scattering phase 
functions. 
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• We show that observations of phase curves of exo- 
planets in reflected hght can be used to constrain 
the underlying scattering mechanisms (e.g. Lam- 
bert versus Raylcigh, etc.) in their atmospheres. 

• We show that observed geometric albedos of cx- 
oplanetary atmospheres can be used to constrain 
their average single scattering albedos, which are 
indicative of their chemical and thermal properties, 
assuming different scattering mechanisms. 

• We demonstrate, using the vector Rayleigh phase 
matrix, how polarization curves can be used to con- 
strain orbital parameters (inclinations and eccen- 
tricities) of exoplanets. 

The simplification of an analytic approach does intro- 
duce some natural caveats. While our requirement of 
semi-infinite atmospheres is reasonable for all gaseous 
atmospheres, from super-Earths to giant planets, our as- 
sumption of homogeneity is meant to represent only an 
average scattering albedo in the atmosphere in the ob- 
served spectral/photometric bandpass. In general, the 
scattering albedo in a gaseous atmosphere would not re- 
main constant at all altitudes and wavelengths. Our ap- 
proach is aimed at providing a simple and efficient proce- 
dure to interpret reflected light observations of exoplan- 



ets which are currently limited by sparse data. 

Detailed numerical models exist in the literature which 
solve the full radiative transfer problem with an inhomo- 
geneous atmosphere in plane-parallel atmospheres. How- 
ever, given the large number of free parameters and com- 
putational complexity of such models, it is generally not 
feasible to use them for formally fitting and retrieving at- 
mospheric parameters from the limited data in reflected 
light. The analytic tools provided in our work makes it 
feasible to derive averaged atmospheric properties, such 
as an average scattering albedo or a representative scat- 
tering phase function, from the data. Our code can be 
used in conjunction with standard parameter estimation 
methods to retrieve average scattering albedos from ob- 
served geometric albedos, as well as orbital parameters 
from observed phase curves and polarization curves of ex- 
trasolar planets. Our code is freely available on request 
by email to either of the authors. 
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